home *** CD-ROM | disk | FTP | other *** search
/ Celestin Apprentice 5 / Apprentice-Release5.iso / Source Code / Libraries / LinAlg 3.1 / LinAlg / fminbr.cc < prev    next >
Encoding:
Text File  |  1995-12-04  |  5.6 KB  |  164 lines  |  [TEXT/ttxt]
  1. // This may look like C code, but it is really -*- C++ -*-
  2. /*
  3.  ************************************************************************
  4.  *
  5.  *              Numerical Math Package
  6.  *
  7.  *            Brent's one-dimensional minimizer 
  8.  *
  9.  *         finds a local minimum of a single argument function
  10.  *              over the given range
  11.  *
  12.  * Input
  13.  *    double fminbr(ax,bx,f,tol)
  14.  *    const double ax                a and b, a < b, specify the interval
  15.  *    const double bx          the minimum is to be sought in
  16.  *    double (*f)(const double x)    Ptr to the function under investigation
  17.  *    const double tol        Acceptable tolerance for the minimum
  18.  *                    location. It is an optional parameter
  19.  *                    with default value EPSILON
  20.  *
  21.  * Output
  22.  *    Fminbr returns an estimate to the minimum location with accuracy
  23.  *    3*SQRT_EPSILON*abs(x) + tol.
  24.  *    The procedure always determines a local minimum, which coincides with
  25.  *    the global one if and only if the function under investigation is
  26.  *    unimodular.
  27.  *    If a function being examined possesses no local minimum within
  28.  *    the given range, Fminbr returns either the left or the right end
  29.  *    point of the interval, wherever the function value is smaller.
  30.  *
  31.  * Algorithm
  32.  *    G.Forsythe, M.Malcolm, C.Moler, Computer methods for mathematical
  33.  *    computations. M., Mir, 1980, p.202 of the Russian edition
  34.  *
  35.  * The function makes use of the "gold section" procedure combined with
  36.  * the parabolic interpolation.
  37.  * At every step program operates three abscissae - x,v, and w.
  38.  *     x - the last and the best approximation to the minimum location,
  39.  *        i.e. f(x) <= f(a) or/and f(x) <= f(b)
  40.  *         (if the function f has a local minimum in (a,b), then both
  41.  *           conditions are met after one or two steps).
  42.  *    v,w are previous approximations to the minimum location. They may
  43.  *    coincide with a, b, or x (although the algorithm tries to make all
  44.  *    u, v, and w distinct). 
  45.  * Points x, v, and w are used to construct an interpolating parabola,
  46.  * whose minimum is to be regarded as a new approximation to the minimum
  47.  * location if the former falls within [a,b] and reduces the minimum 
  48.  * encompassing interval to a larger extent than the gold section procedure.
  49.  * When f(x) has a second derivative positive at the point of minimum
  50.  * (not coinciding with a or b) the procedure converges superlinearly
  51.  * at a rate of about 1.324
  52.  *
  53.  ************************************************************************
  54.  */
  55.  
  56. #include "math_num.h"
  57.  
  58.  
  59. double fminbr(                // An estimate to the min location
  60.     const double ax,        // Specify the interval the minimum
  61.     const double bx,        // to be sought in
  62.     double (*f)(const double x),    // Function under investigation
  63.     const double tol)        // Acceptable tolerance
  64. {
  65.   double a = ax, b = bx;        // Current interval
  66.   double x,v,w;                // Abscissae, descr. see above
  67.   double fx;                // f(x)
  68.   double fv;                // f(v)
  69.   double fw;                // f(w)
  70.   const double r = (3-sqrt(5.0))/2;    // Golden section ratio
  71.  
  72.   assure( tol > 0, "Tolerance must be positive");
  73.   assure( b > a, 
  74.      "Left end point of the interval should be strictly less than the "
  75.      "right one" );
  76.  
  77.   v = a + r*(b-a);  fv = (*f)(v);       // First step - always gold section
  78.   x = v;  w = v;
  79.   fx=fv;  fw=fv;
  80.  
  81.   for(;;)        // Main iteration loop
  82.   {
  83.     const double range = b-a;        // Interval over which the minimum
  84.                     // is sought for
  85.     const double midpoint = (a+b)/2;
  86.     const double tol_act =        // Actual tolerance
  87.         SQRT_EPSILON*fabs(x) + tol/3;
  88.  
  89.        
  90.  
  91.     if( fabs(x-midpoint) + range/2 <= 2*tol_act )
  92.       return x;                // Acceptable approximation is found
  93.  
  94.                     // Compute a new step with the gold
  95.                     // section
  96.     double new_step = r * ( x < midpoint ? b-x : a-x );
  97.  
  98.  
  99.                 // Decide on the interpolation  
  100.     if( fabs(x-w) >= tol_act  )        // If x and w are distinct
  101.     {                    // interpolatiom may be tried
  102.       register double p;         // Interpolation step is calcula-
  103.       register double q;                  // ted as p/q; division operation
  104.                                         // is delayed until last moment
  105.       register double t;
  106.  
  107.       t = (x-w) * (fx-fv);
  108.       q = (x-v) * (fx-fw);
  109.       p = (x-v)*q - (x-w)*t;
  110.       q = 2*(q-t);
  111.  
  112.       if( q > 0 )            // Formulas above computed new_step
  113.     p = -p;                // = p/q with wrong sign (on purpose).
  114.       else                // Correct this, but in such a way so
  115.     q = -q;                // that q would be positive
  116.  
  117.       if( fabs(p) < fabs(new_step*q) &&    // If x+p/q falls in [a,b] and is not
  118.      p > q*(a-x+2*tol_act) &&    // too close to a and b, and isn't
  119.      p < q*(b-x-2*tol_act)  )       // too large, it is accepted
  120.        new_step = p/q;
  121.                     // If p/q is too large then the
  122.                     // gold section procedure would
  123.                     // reduce [a,b] to larger extent
  124.     }
  125.  
  126.     if( fabs(new_step) < tol_act )    // Adjust the step to be not less
  127.       if( new_step > 0 )        // than tolerance
  128.     new_step = tol_act;
  129.       else
  130.     new_step = -tol_act;
  131.  
  132.                 // Obtain the next approximation to min
  133.                     // and reduce the encompassing interval
  134.     register double t = x + new_step;    // Tentative point for the min
  135.     register double ft = (*f)(t);
  136.     if( ft <= fx )
  137.     {                                     // t is a better approximation
  138.       if( t < x )            // Reduce the interval so that
  139.     b = x;                            // t would fall within it
  140.       else
  141.     a = x;
  142.       
  143.       v = w;  w = x;  x = t;        // Assign the best approx to x
  144.       fv=fw;  fw=fx;  fx=ft;
  145.     }
  146.     else                                  // x remains the better approx
  147.     {
  148.       if( t < x )            // Reduce the interval encompassing x
  149.     a = t;                   
  150.       else
  151.     b = t;
  152.       
  153.       if( ft <= fw || w==x )
  154.       {
  155.     v = w;  w = t;
  156.     fv=fw;  fw=ft;
  157.       }
  158.       else if( ft<=fv || v==x || v==w )
  159.     v = t, fv = ft;
  160.     }
  161.   }        // ===== End of loop =====
  162.  
  163. }
  164.